library(tidyverse)
library(haven)
library(labelled)
library(scales)
library(purrr)
library(tidyr)
library(broom)
rm(list=ls())


df<- read_rds("clean-GSS-institutional-trust.RDS")

## make it longer
df_long <- df %>%
  pivot_longer(consci:conlegis,names_to = "institution", values_to = "value")%>%
  #drop_na() %>%
  mutate(institution = as.character(recode(institution, 
                                           "consci" = "Scientific Community",
                                           "coneduc" = "Education",
                                           "conpress" = "Press",
                                           "contv"="Television",
                                           "conjudge" = "U.S. Supreme Court",
                                           "conmedic" = "Medicine",
                                           "conarmy" = "Military",
                                           "conclerg" = "Organized Religion",
                                           "conbus" = "Major Companies",
                                           "conlabor" = "Organized Labor",
                                           "confinan" = "Banks and Financial Institutions",
                                           "conlegis" = "Congress"))) 

means <- df_long %>%
  mutate(value = as.numeric(recode(value, 
                                   "1" = "3",
                                   "2" = "2",
                                   "3" = "1"))) %>%
  group_by(year, institution) %>%
  summarise(mean_trust_i_per_year = weighted.mean(value, wtssps, na.rm=TRUE),
            sd = sd(value, na.rm=TRUE),
            se = sd / sqrt(n()),
            n_non_na = sum(!is.na(value)))

#get total n
sum(means$n_non_na)

ggplot(means, aes(x=year, y=mean_trust_i_per_year)) + geom_point() +
  geom_line(linetype="dashed") +
  theme_bw() +
  facet_wrap(vars(institution)) +
  geom_errorbar(aes(ymin = (mean_trust_i_per_year -1.96*se), 
                    ymax = (mean_trust_i_per_year + 1.96*se))) +
  labs(x="", y="Category",
       subtitle="As far as the people running these institutions are concerned, would you 
say you have a great deal of confidence, only some confidence, or hardly
any confidence at all in them? N=592,648 (GSS)
Mean and 95% Confidence Intervals") +
  scale_y_continuous(limits = c(1,3), breaks=c(1,2,3),
                     labels = c("Hardly Any", "Only Some", "A Great Deal"))
ggsave("PNAS/3.png", width=12, height=8)



###############################################
#diff means test 2000 vs 2004
test <- df_long %>%
  subset(year == 2000 | year == 2024) %>%
  select(value,institution, year, wtssps) %>%
  mutate(value=wtssps*value)%>%
  unique()


mean_values <- test %>%
  group_by(institution, year) %>%
  summarise(mean_value = mean(value, na.rm = TRUE), .groups = "drop")


t_test_results <- test %>%
  group_by(institution) %>%
  nest() %>%
  mutate(
    # Count N for each year
    Ns = map(data, ~ .x %>%
               group_by(year) %>%
               summarise(n = n()) %>%
               pivot_wider(names_from = year, values_from = n, names_prefix = "n_")
    ),
    # Run t-test (unpaired by default)
    t_test = map(data, ~ t.test(value ~ factor(year), data = .x)),
    tidy_result = map(t_test, tidy)
  ) %>%
  unnest(tidy_result) %>%
  unnest(Ns) %>%
  select(institution, estimate, 
         mean24=estimate1, mean00=estimate2, tstat=statistic, 
         p.value,
         conf.low, conf.high, n_2000, n_2024)